% MOC DEMO TO DO:
%  1) find a way to approximately verify the code (w/ MCNP or SN)
%  2) precomputation of various factors
%  3) add a cartesian mesh ability for creating regions (w/o circles for
%     now)  This will be very nice for benchmarking, debugging, and
%     visualizing flux maps
%  4) add reflective boundary condition (cyclic tracking verify???)
%  5) add ability to visualize outgoing angular flux distributions at a
%     global edge
%  6) volume-corrected tracks

%  DEBUG ISSUES
%   * there is a pretty drastic flux change in the same volume cells for
%   difference size cells at the periphery (i.e. 0 1 2 .. 10  gives one
%   central phi while 0 1.5 2 3 .. 9.5 10 gives an answer far different
%   than what the volume changes should give rise to)  This is probably an
%   issue in tracker
%   * for identical problems using MOC and SN, my MOC gives overestimates
%   at the boundary (I assume the SN is right, as it pretty much converges
%   with higher order angular quadratures; to be verified with MCNP)
%  
%


% main file
% solves for the scalar flux in a 2 region pin cell
% based on code from a. hebert
%clear; 

spaces = [50 100 150 200 250 300];
angles = [50 100 150 200 250 300];
polars = [1 2 3 16];

% define the problem
MeshX           = [0 1 2 3 4 5 6 7 8 9 10];
MeshY           = [0 1 2 3 4 5 6 7 8 9 10];
% MeshY           = [0 1 2 3 4 5 6 7 8 9 10];
% MeshX           = [0 0.25 2 3 4 5 6 7 8 9.5 10];
% MeshX           = [0 1 2 3 4 5]*2;
% MeshY           = [0 1 2 3 4 5]*2;
% MeshY = (0:20)/2;
% MeshX = MeshY;
it = 0;
% for s = 1:5
%     for a = 1:5
%         for p = 1:3
       
NumSpace        = spaces(1);%100;
NumAzimuth      = angles(1);%100;
NumPolar        = polars(3);%1;
% region data.  Inner to Outer
% SigmaT          = [ 0.10 0.01 1.00 ];
% SigmaS          = [ 0.03 0.01 0.50 ];
% Q               = [ 1.00 0.00 0.00 ];
SigmaT = ones( (length(MeshX)-1)*(length(MeshY)-1),1);
SigmaS = 0.1*SigmaT;
Q = 1+zeros((length(MeshX)-1)*(length(MeshY)-1),1); 
Q(1)=1; Q(2)=1;
BC              = 0;    % 0=vacuum, 1=reflective (requires cyclic tracking)
MaxIt           = 30;   % Maximum Inner Iterations
EpsPhi          = 1e-7; % relative scalar flux tolerance
% compute the tracks
% figure(2)
track   = trackerCart(MeshX, MeshY, NumSpace, NumAzimuth, NumPolar);
for i = 1:10
    plotterCart(track, MeshX, MeshY, i)
    pause(0.5)
end
% solve or the region fluxes
flux    = mocsolver( track, BC, Q, SigmaT, SigmaS, MaxIt, EpsPhi );
% 
% it = it+1;
% dflux(it,:) = diag( reshape(flux,10,10) );
% 
%         end
%     end
% end

% MOC, 200 x 200, c = 0.1
%     0.7638    0.8857    0.9045    0.9089    0.9100    0.9100    0.9089    0.9045    0.8857    0.7638
%     0.8858    1.0386    1.0637    1.0697    1.0711    1.0711    1.0697    1.0637    1.0386    0.8858
%     0.9044    1.0637    1.0913    1.0980    1.0997    1.0997    1.0980    1.0913    1.0637    0.9044
%     0.9089    1.0697    1.0980    1.1051    1.1069    1.1069    1.1051    1.0980    1.0697    0.9089
%     0.9101    1.0712    1.0997    1.1069    1.1087    1.1087    1.1069    1.0997    1.0712    0.9101
%     0.9101    1.0712    1.0997    1.1069    1.1087    1.1087    1.1069    1.0997    1.0712    0.9101
%     0.9089    1.0697    1.0980    1.1051    1.1069    1.1069    1.1051    1.0980    1.0697    0.9089
%     0.9044    1.0637    1.0913    1.0980    1.0997    1.0997    1.0980    1.0913    1.0637    0.9044
%     0.8858    1.0386    1.0637    1.0697    1.0711    1.0711    1.0697    1.0637    1.0386    0.8858
%     0.7638    0.8857    0.9045    0.9089    0.9100    0.9100    0.9089    0.9045    0.8857    0.7638
% S80, 200x200 fine mesh, c = 0.1
%     0.7134    0.8474    0.8708    0.8767    0.8782    0.8782    0.8767    0.8708    0.8474    0.7134
%     0.8474    1.0182    1.0497    1.0576    1.0596    1.0596    1.0576    1.0497    1.0182    0.8474
%     0.8708    1.0497    1.0845    1.0934    1.0957    1.0957    1.0934    1.0845    1.0497    0.8708
%     0.8767    1.0576    1.0934    1.1028    1.1053    1.1053    1.1028    1.0934    1.0576    0.8767
%     0.8782    1.0596    1.0957    1.1053    1.1078    1.1078    1.1053    1.0957    1.0596    0.8782
%     0.8782    1.0596    1.0957    1.1053    1.1078    1.1078    1.1053    1.0957    1.0596    0.8782
%     0.8767    1.0576    1.0934    1.1028    1.1053    1.1053    1.1028    1.0934    1.0576    0.8767
%     0.8708    1.0497    1.0845    1.0934    1.0957    1.0957    1.0934    1.0845    1.0497    0.8708
%     0.8474    1.0182    1.0497    1.0576    1.0596    1.0596    1.0576    1.0497    1.0182    0.8474
%     0.7134    0.8474    0.8708    0.8767    0.8782    0.8782    0.8767    0.8708    0.8474    0.7134
% error w/r to MOC (%)
%     6.5895    4.3262    3.7283    3.5423    3.4917    3.4917    3.5423    3.7283    4.3262    6.5895
%     4.3275    1.9644    1.3165    1.1288    1.0756    1.0756    1.1288    1.3165    1.9644    4.3275
%     3.7095    1.3154    0.6235    0.4193    0.3628    0.3628    0.4193    0.6235    1.3154    3.7095
%     3.5466    1.1281    0.4193    0.2067    0.1473    0.1473    0.2067    0.4193    1.1281    3.5466
%     3.5003    1.0782    0.3632    0.1473    0.0866    0.0866    0.1473    0.3632    1.0782    3.5003
%     3.5003    1.0782    0.3632    0.1473    0.0866    0.0866    0.1473    0.3632    1.0782    3.5003
%     3.5466    1.1281    0.4193    0.2067    0.1473    0.1473    0.2067    0.4193    1.1281    3.5466
%     3.7095    1.3154    0.6235    0.4193    0.3628    0.3628    0.4193    0.6235    1.3154    3.7095
%     4.3275    1.9644    1.3165    1.1288    1.0756    1.0756    1.1288    1.3165    1.9644    4.3275
%     6.5895    4.3262    3.7283    3.5423    3.4917    3.4917    3.5423    3.7283    4.3262    6.5895
% 
